rm(list=ls())

load("TSCS_data.RData")

# Table 2

lm_antipluralist <- plm(antipluralism_gov_seat_share ~ 
                          lag(antipluralism_gov_seat_share,1:2)  + 
                          citizen_support_lag + 
                          lag(polarization,1) +
                          lag(clientelism_gov_seat_share,1) +
                          lag(region_citizen_support,1) + 
                          lag(region_antipluralism_gov_seat_share,1) + 
                          lag(region_liberal_democracy,1) + 
                          lag(democratic_stock,1) + 
                          lag(GDP_capita,1) + 
                          lag(life_expectancy,1),
                        data = panel,
                        model = "within",
                        effect = "twoways")

summary(lm_antipluralist)

cmodel_antipluralist <- pgmm(antipluralism_gov_seat_share ~ 
                               lag(antipluralism_gov_seat_share,1:2)  + 
                               citizen_support_lag + 
                               lag(polarization,1) +
                               lag(clientelism_gov_seat_share,1) +
                               lag(region_citizen_support,1) + 
                               lag(region_antipluralism_gov_seat_share,1) + 
                               lag(region_liberal_democracy,1) + 
                               lag(democratic_stock,1) + 
                               lag(GDP_capita,1) + 
                               lag(life_expectancy,1) |
                               plm::lag(antipluralism_gov_seat_share, 3:5),
                             panel,
                             effect = "individual",
                             model = "onestep",
                             transformation = "ld",
                             indexes = c("country", "year")
)

summary(cmodel_antipluralist)

# Table model

modslin <- list(lm_antipluralist)
vcm_lm <- lapply(modslin, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))
pval_lm <- lapply(1:length(modslin), function(i) {
  pt(abs(coef(modslin[[i]]) / rse_lm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_wool <- sapply(modslin, function(x) round(pbgtest(x, order=1)$p.value, 3))[c(1,0)] # Wooldridge test for serial correlation

modsmom <- list(cmodel_antipluralist)
vcm_gmm <- lapply(modsmom, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))
pval_gmm <- lapply(1:length(modsmom), function(i) {
  pt(abs(coef(modsmom[[i]]) / rse_gmm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_AB <- sapply(modsmom, function(x) round(mtest(x, order=2, vcov=plm::vcovHC(x, cluster="group"))[["p.value"]][[1]] ,3))[c(1,0)] # AR(2) p-val
model1_H <- sapply(modsmom, function(x) round(sargan(x, weights="twosteps")[["p.value"]][[1]], 3))[c(1,0)] # Hansen test p-value

stargazer(modslin,modsmom,
          covariate.labels=c("Anti-pl. government t-1","Anti-pl. government t-2",
                             "Citizen support t-1","Polarization t-1","Clientelism t-1",
                             "Region citizen support t-1","Region anti-pl. government t-1",
                             "Region democracy t-1","Democratic stock t-1","log(GDP t-1)","Life expectancy t-1"),
          type="latex",
          omit.stat	= c("f","n"),
          se = list(rse_lm[[1]],rse_gmm[[1]]),
          p = list(pval_lm[[1]],pval_gmm[[1]]),
          dep.var.caption	= "",
          dep.var.labels	= "Anti-pl. government",
          add.lines=list(c("N Observations", "1895","3576"),
                         c("N Countries", "99", "100"),
                         c("Wooldridge test (p-value)", model1_wool,""),
                         c("Arellano-Bond test (p-value) ", "",model1_AB),
                         c("Hansen test (p-value)", "",model1_H)
          ),
          no.space = T,
          star.cutoffs = c(0.1,0.05,0.01),
          out = "tab_2.tex"
)
